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Abstract 

An adaptive algorithm, based on residual type a posteriori indicators of errors measured in 
L°°(L 2 ) and L 2 {L 2 ) norms, for a numerical scheme consisting of implicit Euler method in time and 
discontinuous Galerkin method in space for linear parabolic fourth order problems is presented. 
The a posteriori analysis is performed for convex domains in two and three space dimensions for 
local spatial polynomial degrees r > 2. The a posteriori estimates are then used within an adaptive 
algorithm, highlighting their relevance in practical computations, which results into substantial 
reduction of computational effort. 

1 Introduction 

Fourth order parabolic equations and corresponding initial-boundary value problems appear in the 
modelling in areas as diverse as biology, phase-field modelling and image processing to name a few. In 
most cases of practical interest one has to resort to numerical methods for their solution, due to complex 
geometry and/or the presence of non-linearities. 

During the last five decades, finite element methods (FEMs) have been widely used to numerically 
solve fourth order elliptic or parabolic problems; see, e.g., [H H2 EH HE1 El 133] and the references 
therein for earlier works. There are, generally speaking, three families of FEMs developed for fourth 
order problems: conforming, mixed and non-conforming. The classical conforming methods (see, e.g., 
[L3] and the references therein) require the construction of complicated elements with a number of 
degrees of freedom devoted to ensuring (^-continuity across the element interfaces. This results into 
limitations in the applicability of conforming methods on general, possibly irregular, meshes |36) and 
their non-trivial extensions to dimensions three (or higher). Mixed methods (see, e.g., (T2J [16] and 
the references therein), whereby the fourth order operator is first transformed into a system of second 
order operators are widely used in practice, but they require very careful treatment in the imposition 
of essential and natural boundary conditions. Non-conforming methods for fourth order problems were 
first presented by [I] and then further developed in [TH] [TO] 1331 ED] and other works. The key idea in 
non-conforming methods is the use of penalties to ensure convergence into the natural energy space of 
the variational problem, despite finite element basis functions being either just continuous (C°-interior 
penalty procedures; see, e.g., jTB] HO]) or completely discontinuous (discontinuous Galerkin interior 
penalty procedures; see, e.g., 1551 1201121] ) . 

Adaptive FEMs based on a posteriori error estimates has been an active field of research in recent 
years, especially for second order elliptic and parabolic problems. For the case of fourth order elliptic 
problems a posteriori error estimators and indicators have been developed, e.g., in [HI [38] HJ EH [6] 
\13\ |27l EI] . A posteriori bounds and adaptive algorithms for parabolic fourth order problems are 
far less developed in the literature. For instance, the development of adaptive algorithms based on 
various types of a posteriori indicators for the Cahn-Hilliard fourth order parabolic problem can be 
found in [3U EH H3- Error control for variational methods for fourth order parabolic equations has 
been predominantly focused to space-discrete mixed or conforming formulations. The recent work |31) 
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deals with goal-oriented error estimation for the fully discrete Cahn-Hilliard problem. Therefore, the 
development of adaptive algorithms based on a posteriori estimators for fully discrete methods for fourth 
order parabolic problems is still largely an unexplored area. 

Advances in a posteriori error analysis of fully discrete schemes with non-conforming spatial dis- 
cretizations of second order parabolic problems have been recently presented [TH1 [2~5] . In , an adap- 
tive algorithm based on the derived a posteriori estimates is also considered. Local residual a posteriori 
error bounds for semi-discrete conforming and mixed spatial discretizations ffor the Cahn-Hilliard prob- 
lem and the Hele-Shaw flow are presented in [22] • Finally, a posteriori error estimates in an L 2 (H 2 )-type 
norm and adaptive algorithms for fully discrete schemes with discontinuous Galerkin methods for fourth 
order problems are proposed in |39| . The derivation of reliability bounds in [39] is based on the elliptic 
reconstruction framework of Makridakis and Nochetto [32~] ; we also refer to [2"9l [2~5] for some relevant 
extensions. 

This work is concerned with the derivation of a posteriori error estimates in weaker than L 2 (H 2 )- 
norms and their use within an adaptive algorithm for a class of discontinuous Galerkin interior penalty 
methods for a fully discrete approximation of the problem: 

u t + A 2 u = f in Q x (0,T], (1) 
u = Vu • n = in ffix (0, T] and (2) 
u = uq ia Six {0} (3) 

with O C l d , d = 2,3 a convex polygonal domain with boundary dQ. More specifically, we derive a 
posteriori error estimators for the error measured in L°°(L 2 ) and L 2 (L 2 ) norms for a numerical scheme 
consisting of discontinuous Galerkin method in space and simple backward-Euler time-stepping for the 
problem - . The a posteriori analysis is performed for convex domains (as is usual for these norms) 
in two and three space dimensions for local spatial polynomial degrees r > 2. To enable the optimality 
of the a posteriori estimators in the L°°(L 2 ) and L 2 (L 2 ) norms, the elliptic reconstruction framework 
is employed. Moreover, the L 2 (L 2 )-norm analysis employs a special test function construction inspired 
from the a priori analysis of FEMs for wave problems in [2J. Somewhat surprisingly, the use of this 
special testing, in conjunction with the elliptic reconstruction, results into the derivation of L 2 (L 2 )-norm 
a posteriori estimators via a standard energy argument. The efficiency of the a posteriori estimators 
is assessed numerically. The reliability bounds are used within two variants of a space-time adaptive 
algorithm. The adaptive algorithm is able to achieve the same error reduction with far fewer degrees 
of freedom compared to uniform meshes, thereby highlighting the relevance of the derived a posteriori 
estimates in practical computations. The simple model problem - ([3]) appears to be sufficient in 
highlighting some of the challenges in the error estimation and adaptivity of finite element methods for 
more complex fourth order parabolic problems. It appears that the derived a posteriori bounds and the 
respective adaptive algorithms can be modified in a straightforward fashion to include the original dG 
method of Baker [J] and C°-interior penalty methods jTSJ HO] ■ 

The remaining of this work is organised as follows. In Section [2] notation is introduced and some 
standard results needed in the subsequent analysis are recalled. The discontinuous Galerkin (dG) 
method for the biharmonic problem, along with the derivation of posteriori error bounds for the dG 
approximation of the biharmonic problem in L 2 -norm are derived in Section [3| The respective fully 
discrete scheme for the parabolic model problem - is given in Section [4j while Section [5] contains 
the derivation of residual type a posteriori estimates of errors in L°°(L 2 ) and L 2 (L 2 ) norms for the 
fully discrete scheme. The efficiency and reliability of the a posteriori estimators is tested on a range of 
uniform meshes in Section [6] The adaptive algorithm utilizing the a posteriori estimates in a series of 
numerical experiments are also presented in Section [6] Some concluding remarks regarding the results 
and possible extensions are given in Section [7] 

2 Notation and preliminaries 

The standard Hilbertian Lebesgue space is denoted by L 2 (ui), for a domain to C R d , (d = 2,3), with 
corresponding inner product (•, - ,) u and norm || • when uj = O, we shall drop the subscript writing 
(-,-,) and || • ||, respectively. We also denote by H s (to), the standard Hilbertian Sobolev space of index 
s > of real- valued functions defined on to C K d , along with the corresponding norm and seminorm 



|| ■ || S)W and | ■ \ SiUJ , respectively. For 1 < p < +00, we also define the spaces L p (0,T,H s (uj)), consisting 
of all measurable functions v : [0,T] —> H s (l>j), for which 

\\v\\lp(0,T;H'>(ui)) 
IMU~(0,T;X) 

Let T be a subdivision of Q into disjoint elements kgT. The subdivision T is assumed to be shape- 
regular (see, e.g., p. 124 in [15) ) and is constructed via smooth mappings F K : k — » n with uniformly 
bounded Jacobian throughout the mesh family considered, where k is the reference element. The above 
mappings are assumed to be constructed so as to ensure fl = U ki =-t-k and that the elemental edges are 
straight segments (i.e., lines or planes). Note that we also use the expression edge to mean side when 
d = 3. 

The broken Laplacian, A/jU, is defined element-wise by (A/j?/)| K := A(u| K ) for all k G T- 

For a nonnegative integer r, we denote by V r (k), the set of all polynomials of total degree at most 

r, if k is the reference simplex, or of degree at most r in each variable, if k is the reference hypercube. 

We consider the finite element space 

S r := {v e L 2 (n) : v\ K o F K € V r (k), n € T}. (5) 

By r we denote the union of all (d — l)-dimensional element edges associated with the subdivision 
T, including the boundary. Further, we decompose T into two disjoint subsets T — dfl U r; nt , where 

r int : = r\dn. 

For two (generic) elements k + , kT € T sharing an edge e = k + n k~ , we define the outward normal 
unit vectors n + and n~ on e corresponding to dn + and 3k~ , respectively. For functions v : f2 — >■ K and 
q : — > M. d , that may be discontinuous across T, we define the following quantities. For v + :— v\ eC g K +, 
v ~ ■= v\ecdn-, q + := q| e c9K+, and := q\ eCdR -, we set 

{ v } ■= + v ~) : { q } ;= i( q + + q") j [ w ] : = q + n + + q~n~, [q] := q+ • n+ + q • n ; 

if e e 9k n 90, these definitions are modified to {v} := v + , {q} := q + , fuj := v + n, [q] := q + • n. With 
the above definitions, it is easy to verify the identity 

V f vq-nds= [ H-{q}d S + f {v}[q\ds, (6) 
K£T JdK Jr Jr int 

with n denoting the outward normal unit vector on 8k, corresponding to k. 

We define the element size h K := (/^(k)) 1 /^, where /id is the d-dimensional Lebesgue measure; we 
collect the element sizes into the into the element- wise constant function h : — > K, with h| K = h K , 
k E T and h = {h} on T. Also, for two (generic) elements k + , kT sharing an edge e := dn + <~)dK~ C Tint, 
we define h e := fid-i(e)- 

As we shall be dealing with mesh adaptive algorithms below, we assume that all sequences of meshes 
considered in this work are locally quasi- uniform, i.e., there exists constant c > 1, independent of h, 
such that, for any pair of elements k + and k~ in T which share an edge, 

c- 1 < h K +/h K - < c. (7) 

Finally we recall a series of some (standard) results used throughout this work; their proofs can be 

found, e.g., in da nam El HI]. 

Lemma 2.1 (approximation property) Let < m < r + 1 and T be a subdivision of ft, Sic M. d . 
Then there exists a constant C app , independent of h K , such that for any u € H m (tt) and k eT, there 
exists p : C(k) — > E, with p o F K G V r {n) and 

\u - p\ jlK < C app h™- j \u\ m . K , < j < m. (8) 



IKt)||? w dt < +00, for l<p<+oo, 



esssup 0<t<T ||?;(i)|| S:U < +00, for p = +00. 



(4) 



Lemma 2.2 (inverse estimate) There exists a constant Ci nv , independent ofh K , such that 

\p\j,n < C mv h l ~ ] \p\i, K , < i < j < 2, (9) 
for all p : C(k) R, with p o F K £ V r (n). 

Lemma 2.3 (trace inequality) For every u £ H (k), with k £ T, there exists a constant C tr > 
independent of h K such that 

Mlo^Ctrih-^Mln + hMln)- (10) 

Lemma 2.4 (Poincare-Friedrichs inequality There exists a constant C p f, independent ofh K , 

such that for any u £ L 2 (fl), with u\ K £ H 2 (k) for all k £ T, we have 

\Hln + Hln < C pf (\u\l n + ||h- 3 / 2 M||2 T + Wh-^lVu}]^) . (11) 

3 Discontinuous Galerkin method for the biharmonic problem 

We consider the biharmonic equation 

A 2 u = cj) inn, (12) 
with homogeneous essential boundary conditions 

u = , Vu • n = on d£l, (13) 

where n denotes the unit outward normal vector to d£l and (f> £ L 2 (Vl). Then the regularity of the 
problem implies that u £ iJ 4 (fi) n i?o( fi ) H3- 

Upon defining the lifting operator C : S := S r + H 2 (fl) —> S r by 

C{v)ipdx = /"(H-{V^}-{^}[V^ds Vip£S r , (14) 



the (symmetric) interior penalty discontinuous Galerkin (dG) method for (12), (13 I is given by: 

find £ S r such that Biu^Vh) = l( v h) V«/, £ S r , (15) 
where the bilinear form B:cSx5->K and the linear form I : S — > R are given by 



with 



B(w, v) := ^ (A h wA h v + £{w)A h v + A h w£(vjj dx + B p (w, v) (16) 

B p (w,v) := J (jH • [u]+£[Vw][Vw]) ds, 
and 

l(v) := j cj)vdx, (17) 



n 



respectively, for w,v £ S. The piecewise constant discontinuity penalization parameters er, £ : T — > R 
are given by 

a\ e = ^(hle)' 3 , e|e=^ (h|e)- 1 , (18) 



respectively, where ao > and £o > 0. To guarantee the stability of the IPDG method defined in (15), 
Co and £o must be selected sufficiently large. 

Note that this formulation is inconsistent for trial and test functions belonging to the solution space 



S. However, when w,v £ S r , in view of (14), (16) gives 



B(w,v)= / A h wA h vdx+ / ({VAto} • [«] + {VAw} • {wj 

Jn Jr v (19) 

- {Aw}[Vv] - {Av}[Vw} +<t[w] ■ H +C[Vu;][Vw]) ds; 



therefore, (15) coincides with the symmetric version interior penalty method presented in [37] . For the 
bilinear form £?(•, •) in ( 16 1 we have the continuity and coercivity with respect to the energy norm on iS 
defined by 

HMD = (HAHIn + II^MIIr + llVltVHIIr) 1 - ( 20 ) 



Lemma 3.1 ([23j) For sufficiently large do > and £o > there exist ■positive constants C con t and 
C coer , depending only on the mesh parameters such that 



\B(u,v)\ < C co „t|||u||| lll^lll Vu,v € S and 
B(u,u) > C coer |||w||| 2 Vm G S . 



(21) 
(22) 



An a posteriori bound for the energy norm error of the dG method (15) for (12) - (13) has been 
considered in [24] . Now, we shall present an a posteriori bound for the L^-norm error (cf. [35 for a 
corresponding result for the second order problem). 

Theorem 3.2 (L 2 -a posteriori bounds for the elliptic problem) Let u € H A (Vl) n Hq(JI) be the 



solution of | 12)-mM, iih € S r be dG approximation (15) associated with the mesh T ■ Then, there exists 
a positive constant C pT^ , independent ofT, h, u and Uh, such that 

\\u-u h \\<E(T,u h ,(f>), (23) 

where 

£ (T, u h , <j>) := C B (||h 4 ~ A / 2 (0 - A 2 u h )\\ 2 + ||h( 7 - A )/ 2 [VA^]|| 2 in 
+ E i 1 + &) ll[v^]f + h\r x (1 + a 2 ) \\{u h \ 



|h( 5 - A )/ 2 [A%]||^ 



(24) 



and A := 2(2 - irrm{2,r- 1}). 
Proof The dual problem 



A 2 z = u — Uh ='■ e in Q, 



(25) 



with homogeneous essential boundary conditions z = Vu • n = on dfl clearly satisfies e € L 2 (Q) and, 
therefore, the following regularity estimate holds 



<C„ 



(26) 



Using p5j ), integrating by parts twice, applying ^ and (14) as well as the regularity of the dual solution, 
z e H 4 (n), we have 



lell 2 



= A 2 zedx = / A h zA h edx- / [Ve]{Az}d,s + / [e] • {VAzjds. (27) 

k£ j-Jk. Jq. Jr Jr 



By using the fact that u is a weak solution and integrating the term involving AhzAhUh by parts, we 
arrive at, 



lei 2 



B[u,z)- / A h zA h u h dx + / [\7u h ]{Az}ds- / [u h j ■ {VAzjds 
Jn Jr Jr 

l{z)- / zA 2 u h dx+ / {z}[VA h u h ]ds- / {Vz} ■ [A fe fi fc ] ds 
+ y [V%]{A4ds - y [%1 • {VAz}ds. 



(28) 



Using the standard orthogonal L 2 -projection, II : S — >• S r , of z , we can derive the following identity by 
integrating by parts and using ^ and ( 14 1 as follows, 

= l(-Uz) - B(u h , -ILz) 



J2 I ((cb-A 2 u h )(-Uz)-C(u h )A h (-Uz))dx + f {-Uz}[VA h u h ] ds 

- f { V(-ILs)} • {A h u h j ds - f (aluhj ■ [-Hz] + £[Vu h ] [V(-Ilz)}) ds. 
Jr int Jr 



(29) 



Using ( 14 ) in ( 29 ) and combining ( 28 ) and <|29h , we get 



le|| 2 =||e|| 2 



(30) 



l(-Hz) - B{u h ,-H.z) 
= I (<p- A 2 h u h ){z - nz)dx + / ({z- Uz}[VAu h ]ds - [Au h ] ■ {V(z - ELz)})ds 

- JlM ■ ({VA(z - Hz)} + (7 h- 3 [z - Hz]) ds 

+ |[V(u k )] ({A(z - nz)} + eoh-^v^ - n*)])ds. 

The assertion then follows by applying Young's inequality, the trace inequality (10) where appropriate, 
the approximation property ((8l and the regularity of the dual problem on each of the terms on the right 



hand side of ( 30 1 . 



□ 



Remark 3.3 If a smooth C 1 subspace of the finite element space exists, such as Argyris elements in two 
dimensions, or corresponding constructions in three dimensions, it is possible to establish an a posteriori 
L 2 bound without dependence on penalty parameters; indeed, these terms would vanish from \3(fy if the 
projection, II, could be defined onto the smooth subspace of S r . 



Remark 3.4 It is interesting to note that the a posteriori error bound of (23) reflects the suboptimal 
L 2 -norm error convergence of the dG method when quadratic polynomials are applied. Similar behaviour 



is observed theoretically and numerically in \23\/ and in \31^ in the context of the a priori error analysis 
of the same method. 



4 DG method for the parabolic problem 

Throughout the remaining of this work, we shall denote by u the weak solution of the problem ([l])-(|3| 
in variational form: find u G £7^(0, T; i? 4 (fi) n H$(Q)) such that 

{u t ,4>)+B{u,<j>) = {f,<t>) V0eff 2 (fi), 

L 2 (n) in n x {o}. 

We consider a subdivision of the time interval (0, T] to be the family of intervals {(t n ~ 1 ,t n ] ; n = 
1,...,7V, with t° = 0, t n - x < t n and t N = T} , with local time-step A„ := t n - t n ~ r . Associated 
with this time-subdivision, let T n , n — 0,...,N, be a sequence of meshes which are assumed to be 
compatible, in the sense that for any two consecutive meshes T n -i and T n , T n can be obtained from 
Tn—i by locally coarsening some of its elements and then locally refining some (possibly other) elements. 
The finite element space corresponding to T n will be denoted by S r,n and the respective dG bilinear 



form by B n {■,■). The backward Euler-dG method for approximating (31) is then given by: for each 
n = 1, . . . , N, find 

U n e S r - n such that (- ,V) + B n (U n , V) = (/", V) W e S r > n , (32) 

where /°(-) := /(-,0) and /"(•) for n = 1, . . . , N is a piecewise polynomial of degree p in time L 2 - 
projection in time of the source function /. In practice, it suffices to take p = to achieve a first- 
order-in-time convergent method. We also set U° := II uo, with II : L 2 (il) — > S r '° is the orthogonal 
L 2 -projection operator onto the finite element space S r '° . 



5 A posteriori bounds for the parabolic problem 

We shall derive a posteriori error bounds for the backward Euler-dG method |32| ) measured in L°°(L 2 )- 
and L 2 (L 2 )-norms. To this end, we shall employ an energy argument (with carefully defined test 
functions) in conjunction with the elliptic reconstruction technique [32 , 29, 25J. 

We begin by extending the sequence {U n } n=1 N of numerical solutions into a continuous piecewise 
linear function of time 

t - t n ~ x t n - t 

[/(0)=II%) and U(t):= r U n + — U n ~ l (33) 

A,i A„ 

for t £ (t n -i,t n ] and n = 1, . . . , N. Further, the discrete elliptic operator A n : S r ' n — > S r ' n is defined by 

for^eS^, (A n <p,x)=B n (<t>, X ) V X eS r '". (34) 

We now give definitions of the estimators involved in the estimation of the parabolic part of the 
error. Estimators at time step n are denoted by oo, n subscript and 2,n subscript will be used for the 
cases of L°°(L 2 )- and L 2 (L 2 )-bounds presented below, respectively. 

Definition 5.1 (estimators for the parabolic error) We define the coarsening or mesh-change es- 
timators by 

1 n— 1 

7oo, n := ^-|| (/ -rnir 1 - 1 !! 2 , 72 , n : = ||(/-n")f/"- 1 || 2 + ^||(ff-n i; - 1 )c/ J - 1 || 2 , (35) 

n i—l 

the time-error evolution estimators by 

n-l 

Voo, n ■= \\9 n -9 n - l \\ 2 K, m, n ■■= \\g n -g n ^ 1 K + Y l ^y-9 i - 1 \\ 2 , (36) 

i=l 

g n := A n U n — n n /" + /"; the data approximation error in time estimators by 

A»,»:= f ||/ n -/|| 2 di, p 2 ,n~K[ t \\f l ~f\\ 2 dt, (37) 

and an additional space estimator given by 

»*»,« ~£(f n ,U n -U n - 1 ,g n -g n - i y, (38) 

where T n := T n n T n -i finest common coarsening of T n and T n -i for each n = 1, . . . , N . 
Using the notation above, we are ready to state the main result. 



Theorem 5.2 (a posteriori bound) Let u £ L 2 (0, T; H 4 (£l) n H$(Q,)) be the solution of pip , U 
be the approximation obtained by the dG method (32) and defined by (33). Then there exist positive 
constants and Ci, independent ofl~ n , h, u and U , for any n = 1, .. . , N such that 



N 



< Coo M|e(0)|| + ,71 H~ Voo.n "h fioo.n a I 

+ fVf? 0o , n ] + max {E(T n ,U n ,g n )}) (39) 

IMU 2 (0,T;L2(fi)) < C 2 ['||e(0)|| + (7 2 ,n +J/2,„ + /?2,n) A„^ 

+ [J2 £(T n ,U n ,g n f A„] ]. (40) 



\n=l 

/ W 



The proof of this theorem will be the content of the remaining of this section, split into a number of 
intermediate results. 

We begin by defining the elliptic reconstruction uj n £ Hq(CI), of U n to be the solution of the elliptic 
problem 

B n (cj n ,v) = (g n ,v) yveH^n) (41) 

where, as above, g n :— A n U n — II n / n + We note that under the assumptions on the domain 
0, we also have w™ g H 4 (il). We also extend the elliptic reconstruction into a continuous piecewise 
linear-in-time function 



w(i) := 



t-t n - x „ t n -t 



-LO 



(42) 



An A n 

for t e (t n —i, t n ] and n = 1, . . . , N. Finally, we introduce the error decomposition 

e := U — u = p — e, where e :— lo — U, and p := ui — u, (43) 
where p and e are understood as the parabolic and elliptic error, respectively, and we set e n := e(t n ). 
Lemma 5.3 (Error identity) For all t £ (t n -i)t n ], n = 1, 2, . . . , N, we have 



(p t ,v)+B n (p,v) = (e u v) + ((I-U n )U t ,v) + 



t-f 
A,, 



(g n ~g n -\v) +(f»-f,v), 



(44) 



for any v G Hq(CI), with I denoting the identity mapping. 



Proof Firstly, from (41 1 and (34) we have 

B n (w n ,v) - (f n ,v) = B n (u n ,n n v) - (n n /™,n™w). 

Also, using the method (32 1 and the definition of the L 2 -project ion we deduce 

{U t ,v) = ((I-n n )U t ,v) + (U t ,U n v) = ((I-n n )U t ,v)- (B n (U ni U n v) - {U n f n ,Il n v)). 
For the elliptic reconstruction error we also have, 

t-e 



B n (uj-u} n ,v) 



A, 



(g n -g»-\v). 



Lastly, for the terms on the left hand side of (44), we compute 

(e t ,v) + B n (p, v) = {U t ,v) + B n (uj, v) -((«*,»)+ B n (u, v)) = {U u v) + B n (cj, v) - (/, v). 



(45) 



(46) 



(47) 



(48) 



Using (45), (45), and (47) in (48), along with the identity e = p — e completes the proof. 

The a posteriori bounds (39) and (40) will be derived by selecting special test functions v in the 
energy identity ( 44 ) above, along with estimation of the terms on the right-hand side of ( 44 ) . More 
specifically, we consider the following two test functions, v := p for the L°°(L 2 ) case, and 



v(t, 



J P (s,-)ds, te[0,T], 



(49) 



for the L 2 (L 2 ) case; this choice is motivated by Baker [3], who used a similar construction for the proof 
of a priori bounds for the second order wave problem. The latter test function has most notably the 
following properties: 



v E H 4 (fl) n H 2 (n) as p e H 4 {n) n H$(n) a.e. in [0,T], 
v{T, ■) = = Av(T, ■), Vu(T, -) = 0, and 
v t (t,-) = -p{t,-), a.e. in[0,T]. 



(50) 



Next, we consider two auxiliary functions which are needed in the consequent proofs. More specifi- 
cally, on each interval t £ (t n -i, t n ], for n = 1, . . . , N, we define 

G(t) := (I-n n )J7 + ^ n , with ^ n :=-(/-n B )Z7 n - 1 + V n ~ 1 ) *P° ■= 0, (51) 

and 

\ / f — f n \ 2 \ 

n n I >> L \ , ri „n~\\ i a n „,;^i, on ^ n I „n „n-l\ , an-1 qO 



G ®--=Y{-^— ) (9 n -g n - x )+B n , with 9 n :=-^(g n -g n -')+0 n -\ 6" = 0. (52) 

We note that, for each n = 1, . . . , N, we then have G(t n ) = ip n , G(t n ) = 9'\ 

t — t n 

G t (t):=(I-n n )U t , and G t (t) := — — (g n - g n ~ l ) . (53) 

An 



The following estimates will be used in the proof of Theorem 5.2 
Lemma 5.4 Let t 6 (0,T]. Then, we have 



/ <G t ,p)dt < ^iKn"-/)^- 1 !! max T ||p|| 

/•T W 

/ (Gt,p)df < E^l^-S"" 1 !! max 

■ y " n=l 



(54) 



N 

HI (55) 



0<t<T 

(f n -L P )dt < J ||r-/||dt max T ||H|. (56) 

Proof The proofs of these estimates are immediate via Cauchy-Schwarz-in-space and Holder-in-time 
inequalities. 

□ 

In the following three lemmata, we prove bounds for the corresponding terms to the ones in Lemma 



5.4 when testing with v given in (49). 



Lemma 5.5 With the above notation, we have 

N »t" N n-1 i „ t " i 

]T / ((j-n-)t/ t ,z)} < ^ (A„||(j-n-)c/— 1 || 2 +A ri ||^(n i -li^- 1 )^- 1 !! 2 ) 2 ( / || P || 2 d<) 2 . 

7»=1 1 n=l i=l 1 

(57) 

Proof Recalling the definition of G, an integration by parts with respect to time gives 

N t*> N n N t" N i n 

W ((I-n n )U t ,v)dt = ^\(G(t),v(t))]\ + J2 [ (G,-v t )dt = J2 (G,-v t )dt. (58) 



n=l 1 n=l n=l 



We recall the properties of d in (50) and we estimate theright-hand term further: 



E / (G,-v t )dt<J2( ||G|| 2 dt) 2 (/ ||H| 2 dt) 2 . 



(59) 



The assertion then follows by estimation of the time integral of ||G|| 2 : 

||G|| 2 di < A„||(/ - n")^"- 1 !! 2 + A„||^-i|| 2 , (60) 

1 

and noting that ip n _ i = J27=i (H 4 - IP - 1 ) 17*' 1 . 



Lemma 5.6 With the above notation, we have 
N r* n f-+ n N / n_1 x \ * / r* n \ • 



7i— 1 ~ ' -,t n— 1 i—1 

Proof Recalling the definition of G, an integration by parts with respect to time gives 



" pt" , _ ,n " ft" 

E / ^" - 9 n_1 ^) d * - E / di - 



(61) 



(62) 



We recall the properties of v in ( 50 1 and estimate the right-hand side further: 

E/ (G,-v t )dt<J2( \\G\\ 2 dt) (/ ||p|| 2 dt) . (63) 

The assertion then follows by estimation of the integral of ||G|| 2 : 

rt" 

V 3 1 1 „n ^n— 1 II 2 



||G||^<A^|! 3 "- 5 ™- i ||^ + A„||^_ i r (64) 

-i 

and noting that n _i = _ y (5* ~ ff^ 1 )- 

Lemma 5.7 W^it/i i/ie above notation, we have 

E / </" - /> ") dt ^ c «w E ( / A "ii/" - /n 2d< ) ( / iHi 2d *) • ( 65 ) 

Proof As /" is the L 2 -projection of / in time, we have 

E / (/" - /, v)dt = E / (/" -f,v- C)dt, (66) 



for the lowest order time approximation C n (") : — ^n, 1 Itn-i v(t, -)dt of v. With the approximation 
property of ( n in time, we deduce 

\\v-C\\ 2 dt<Cl w Xlf \\v t \\ 2 At, (67) 
and recalling that v t = — p, the Cauchy-Schwarz inequality implies 

E/ (f-mdt<c app j2( \i\\r - f\\ 2 dt) (/ yi 2 dt) . 

To complete a posteriori error bounds in Theorem |5.2[ we also need the following two lemmata in 
which the elliptic error terms e and et are estimated by fully computable residuals. 



Lemma 5.8 Let e be as in (43). Then, we have 



/ Ikll 2 dt < 2 -^£(T N ,U N ,g N ) + E ^£(T n ,U\g n Y 

J 6 n =l 6 



(69) 



Proof Noting that ((f - t n_1 )/A„) 2 < § and ((i n - t)A„) 2 < §, we have 

rT N 



I llefd^^^rile^p + He^xll 2 ). (70) 

J ° n=l 6 V ' 



The assertion then follows by Theorem 3.2 



□ 



Lemma 5.9 Let e &e as m |^<3) and r 6 [0,T]; then, we have 

N 



/ (e t ,P>dt< VfC^.tP-CT'- 1 , 

J ° n=l 



g n - g"- 1 -) waxM, (71) 



where 7~ n '■= T n (1 7^-i denotes the finest common coarsening ofT n and T n —\, n = 1,...,N. 

Proof We have e t (t) = (e„ — e n _i)/A n , for t e (t n _i,t„] and n = 1,...,N. Denoting r := t r+ i/2 and 
r := max{fc : ij. < t, fc = 1, . . . , iV}, we then have 

/ (e t ,p)dt= V / — (e„-e n _i,p)dt< max JIpII V ||e n -e n _i||. (72) 

JO n=1 it"- 1 A n 0<t<T ^ 



We now observe that the finite element function z in the proof of Theorem 3.2 can be selected from 
a subspace of S r : in particular, we can select the finite element subspace corresponding to the finest 
common coarsening mesh 7~ ni for n = 1, . . . , N. Then, following completely analogous argument as in 
the proof of of Theorem |3.2| we can arrive to the bound 



□ 



e»-i|| <£(T n ,U n -U n - 1 ,g n -g n - > 

which already yields the result. 



Remark 5.10 Note that the following simpler alternative bound for the term in Lemma \ 5.9\ is also 
possible, 

/ (e t ,p)dt< £f (T n ,U n ,g n ) max .\\p\\. (73) 



5.9 



to the 



This bound shifts the emphases from the finest common coarsening mesh, 7~ n , in Lemma 
elliptic estimators acting on meshes at each time step only which can be of practical importance when 
implementing adaptive algorithms based on the estimators. 



Proof of Theorem |5.2| To conclude the proof, we estimate the left-hand side of (44 1 in each case of 
the test functions: V to derive L 2 (L 2 )-norm a posteriori bound and p for the i°°(L 2 )-norm bound. First 



we deal with the L (L ) case; we start by integrating (44 1 by parts in time 



{e t ,v)+B(p,v)dt = / {e,-v t )dt+[{e,v)]l - / B(v t ,v)dt 
ii Jo Jo 

T r T r T 



f (p,p)dt-J (e,p)dt-(e(0),v(0))- J \j t B &v)dt (74) 
= jT \\p\\ 2 dt - £(e, p)dt - <e(0), «(0)> + \b{v(0), 0(0)). 

We also have 

<e(O),0(O)) < ||e(O)||||0(O)|| < ||e(O)||C pfJ B(0(O),0(O)). (75) 



Using ( 74 ) and ( 75 ) in ( 44 ) after integration over each interval (t n , t n ] and summation with respect 



to n, we get, 

2 

\L 2 (0,T,L 2 ((l)) 



IIpII 2 



N t « 

<l|e(0)|| 2 + ^ jf^ ((e,p)+((I-IL n )U u : 



t-tf 



{g n -g n -\ v)+(r-f, e) dt. (76) 



The bound ( 40 1 now follows upon using the triangle inequality 

ll e IU 2 (0,7\L 2 (fi)) < IIpIIl 2 (0,T,L 2 (O)) + ll e llL 2 (0,T,L 2 (O)) ; 



Young's inequality and Lemmata 5.5 5.6[ 5.7 and 5.8 



For the L°°(L 2 )-norm case, upon testing with v — p, we deduce for the left-hand side of (44 1 by 
integrating by parts to some r £ [0,T] 



(e u p)+B(p,p)dt=\\p(r)\\ 2 -\\p(0)\\ 



(e*,P>dt- 



B(j>, P )dt. 



(77) 



Choosing r such that ||/j(t)|| = max <t<T ||p||, using the triangle inequality, ||p(0)|| < ||e(0)|| + ||e(0)|| 
and (n\ in (Eil, we get, 



\\p\\U(o,T,mn )) + B( PlP )dt<\\e(0)\r + \nOW+ ((e t ,p) + (G t ,p) + (G t ,p) + (f n -f,p))dt (78) 



where G and G are given by (52) and (51). The bound (39) follows again using triangle inequality 



IN|L°°(0,T,L 2 (r2)) < IIpI|l°°(0,T,L 2 (O)) + ||e||L=°(0,T,L 2 (fi)), 



(79) 



Lemmata 



5.4 



and 



□ 



5.9 as well as maxo<t<T ||e|| < maxo<f<T £ (T n , U n ,g n ). 

We note that a posteriori bounds in the L 2 (H 2 )-norm of the error have been already considered 
in |39j . along with their application within an adaptive algorithm. The L 2 (_ff 2 )-norm theoretical and 
numerical results appear to be of the expected order of convergence; they are omitted here for brevity. 



6 Numerical Experiments 

For t £ [0, 1] and O := (0, l) 2 , we consider two benchmark problems for which uq and / are chosen so 
that the exact solution u of problem (31) coincides with one of the following solutions: 



Ul (x,y,t) =sin(7rf) 10 2 sin 2 (7ra;) sin 2 (7r 2 /) e - 10(a;2+!/2) , 
u 2 (x,y,t) = sin(207rf) sin 2 (7ra) sinV^e-^+A 



(80) 
(81) 



Solutions u\ and ui are both smooth but 112 oscillates much faster where as u\ exhibits greater space 
dependency of the error. They are defined so as to emphasize different aspects of the estimators at 
hand. Similar examples have been studied elsewhere, for example in |39] in the context of L 2 (H 2 ) -norm 
a posteriori estimators; see also [59 l 150"! 125] for similar examples in the context of second order problems. 

For the numerical experiments, the library FEniCS (http://fenicsproject.org/) was used. For 



each of the examples, we compute the solution of (|32j) using quadratic simplicial finite element spaces 
and with interior penalty parameters cr = £0 = 20 in ( |18[ ), which is sufficient to guarantee stability of 
the numerical scheme. The interior penalty parameters have a known effect on the effectivity indices, 

cf., EES Eg. 

We study the asymptotic behavior of the indicators by setting all constants appearing in Theorem 
5.2 equal to one. We monitor the evolution of the values and the experimental order of convergence 
of the estimators and the error as well as of the effectivity index over time on a sequence of uniformly 



refined meshes with h K 



2 -i/a- 



1. 



,5, k £ T with fixed time steps A 
al order of convergena 

positive quantities a(i) defined on a sequence of meshes of size h(i) by 



max K hz, and 



A ~ max K h 2 K . To this end, we define experimental order of convergence (EOG) of a given sequence of 



(82) 



the accumulated coarsening or mesh change estimators by 



E 



coarsen, oo,m 



m 

'■= ( ^2 7oc,«<0 
n=l 



and E 



coarsen, 2, m 



n=l 



n=l 



and E t 



accumulated time error evolution estimators by 

m tci 
■^timc,oo,m ■ — ( ^ ^ (^oo,n /^oo, n)^n 

accumulated space error estimators by 

E spa ce,co, m := max {£ (%, U n ,g n )} and E spaco , 2 , m 

0<n<N 



(83) 



(^(%,n+/32,n)A„) 5 , (84) 



(85) 



and the inverse effectivity index 

IM|L~(0,t m ;L 2 (Q)) 



E 



timc,oo,m 



- E 



or 



||e||L 2 (0,t m ;L 2 (a)) 



space, oo, m 



^timc,2,m 



■E 



space, 2, m 



(86) 



for the case L°°(L 2 ) and L 2 (L 2 ), respectively. The IEI conveys the same information as the (standard) 
effectivity index and has the advantage of relating directly to the constants appearing in Theorem |5.2| 
The results of numerical experiments on uniform meshes, depicted in Figures 1-4, indicate that the 
error estimators are reliable and also efficient which can be seen from the effectivity index behaviour 
and the EOC of the error and the time and space estimators for both L 2 (L 2 )- and i°°(L 2 )-norm a 
posteriori bounds. 

To further evaluate practical aspects of the derived a posteriori estimators, they are incorporated 
within in two adaptive algorithms; these are outlined in pseudocode as follows 



ImplicitTimeStepControl 

Input: U , f, TOLt imc , min , TOL timc , TOL spaco , . . . 

TOL coarse , Ao, to, T, . . . 

To, Zrefine, SpaceAdaptivity, . . . 
InitialSpaceAdaptivity 

{ Initial condition interpolation and mesh refinement } 
(Uo, 7o):=InitialSpaceAdaptivity(l/o, /, 7o, fre/ine)- 
{Initialize.} 

Set: n = X, X n — A n _i. 
While (t„ < T) 



Set: E tin 



TOL tin 



1 



While (E timc > TOL timc ) 

t n . — t n — \ -f- \ n 

Set: Tt ■= T n 

(U n ,Tn) := SpaceAdaptivity(t/ n _i, . . . 

/, TOL space ,TOL 

coarse j * ■ • 
An ,t n ,T, T'n — l , ^refine) 

compute Etimc • 
if (E timc > TOL timc ) then 

{Shorten timestep.} 

A n ; = A n /2 

Set: Tn ■= Tt 
endif 
End While 

An + 1 - = A n * 2 

n := n + 1 
End While 
Output: U n 

where SpaceAdaptivity (and InitialSpaceAdaptivity) are performed using a standard Dorfler 
marking strategy expressed in pseudocode as follows 



ExplicitTimeStepControl 

Input: Uo, f, TOL t i meimin , TOL timc , TOL spacc , . . . 

TOL CO \o,to,T,... 

To, fine, SpaceAdaptivity, . . . 

InitialSpaceAdaptivity 
{ Initial condition interpolation and mesh refinement } 
(£/o,7b):=InitialSpaceAdaptivity(f/o, /, To, define)- 
{Initialize.} 

Set: n = 1, A n = \ n -l and tn = t n -i + A„. 
While (t n < T) 

(Un,Tn) ■= SpaceAdaptivity(l/„_i, . . . 

/, TOL space ,TOL 

co a rsc j • ■ • 
An ,t n ,T, Tn — 1 , ire fine) 

compute E t i mc . 
if (E timc > TOL timc ) then 

A„+i := A n /V2 
elseif (Etimc < TOL t i m c,min) then 

A„+i := A n * v2 
endif 

tn + l ; — tn "f" A n 

n := n + 1 
End While 
Output: U n 



Space Adapt ivity 

Input : U n — \ , f ', TOL s pacc i TOL coarsc , T n ,t n , X", Tn — 1 j ^refine 

Set: 7n != 7n— 1- 

TJi := SpaceCoarsening(C/ n _i,TOL CO arse,T n ,7^) 
{Refinement} 

compute local elliptic estimators, (LocalEst njK ) K g7- re . 

sum up local estimators and set Sum to tal : = Y2 K eT LocalEst n:K , and compute E spaco . 
While (E apace > TOL spacc ) 

sort (LocalEst n|K ) K e7^ in descending order, set Q := 0. 

Set: Sum = 0. 

While ((Sum < def ine * Sumtotal) and (k£T»)) 
{Dorfler marking } 

Sum := Sum + LocalEst niK 

if (Sum < £ re fi ne * Sum total ) 

Mark k for refinement; Q := {k} U Q. 
End While 

Refine all elements in Q to obtain new mesh T n . 
Solve I n U„--L. 

Solve |32l for U n with U n U n -i, 11™/™, r n and t n on T n . 
compute local elliptic estimators, (LocalEst njK ) reg 7- n . 

sum up local estimators and set Sum tota i := 5Z K eT LocalEst n K , and compute E spaco . 
End While 
Output: U„,T n 

The refinement ratio < £re/ine < 1 and the tolerances TOL spacc > 0, TOL spaco > and TOL coarso > 
are predefined quantities. The value of Refine '•= 0.75 was used throughout the experiments in adaptive 
algorithms. Note that the coarsening tolerance, TOL coarso , (as well as the tolerance for the alternative 



space estimator in L°° case of Remark 5.10) had to be determined experimentally for given space and 
time tolerances and depending on an example. 

The results of experiments with adaptive algorithms as well as a comparison between the two algo- 
rithms, are detailed in Figures 6-8 where we monitor time step size, accumulated degrees of freedom and 
error evolution in comparison to the uniform approach leading to the desired tolerance. The results of 
these test cases imply substantial reduction in degrees of freedom by both adaptive algorithms in order 
to reach the same error tolerance as compared with the uniform approach. This implies a potential 
efficiency gain in solving PDE problems addressed in this work. 

The estimators presented here are found to be suitable for both adaptive time stepping algorithms 
due to their good separated scaling properties in time and in space. The numerical results appear to be 
less sensitive to mesh change, compared to the same adaptive algorithms based on the L 2 (ff 2 )-norm a 
posteriori error estimators presented in |39) . For instance, terms involving (<?" — <?™ _1 ) which is sensitive 
to mesh change (coarsening as well as refinement) scale down sufficiently fast with the a posteriori 
estimators presented in this work, resulting to robust error reduction in an adaptive algorithm. 

Finally, we note that the considerably more computationally efficient ExplicitTimeStepControl 
algorithm (due to absence of time step searching step) was found to reach desired error tolerances (even 
though this is not guaranteed in general) in these numerical experiments. 



7 Concluding remarks 

Residual type a posteriori estimates of errors measured in L°°(L 2 )- and L 2 (L 2 )-norms for a numerical 
scheme consisting of implicit Euler method in time and discontinuous Galerkin method of local poly- 
nomial degrees r > 2 in space for linear parabolic fourth order problems in space dimensions 2 and 3 
are presented. Numerical experiments confirming the practical efficiency and reliability of the a poste- 
riori estimators are also presented, along with the use of these a posteriori estimator within adaptive 
algorithms, ft appears that the derived a posteriori bounds and the respective adaptive algorithms can 
be modified in a straightforward fashion to the original dG method of Baker [?] and to the C°-interior 
penalty methods of [T51 [TO] . Moreover, second order operators can be included in the present analysis, 
as was done in [39 . An extension of these results to nonlinear fourth order parabolic equations remains 
a future challenge. 
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Figure 1: Example (80) with solution ui, L 2 estimators and error. 
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(a) Example (80) solution with fixed time step A ~ h . Inverse effectivity index is asymptotically tending towards value 
0.004. 
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(b) Example (80) solution with fixed time step A ~ h . Inverse effectivity index is asymptotically tending towards value 
0.004. 



Figure 2: Example (80) with solution ui, L°° estimators and error. 
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(a) Example (80) solution with fixed time step A ~ h? . Inverse effectivity index is asymptotically tending towards value 
between 0.005 and 0.006. 
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(b) Example (80) solution with fixed time step A«A . Inverse effectivity index is asymptotically tending towards a value 
near 0.005. 



Figure 3: Example (81) with solution 112, L 2 estimators and error. 
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(a) Example (81) solution with fixed time step A ~ h . Inverse effectivity index is asymptotically tending towards value 
0.004. 






1 1 1 r 




0.2 0.4 0.6 0.8 1 

EOC(E t 2 m ) 



0.2 0.4 0.6 0.8 1 

EOC(E s 2 m ) 



0.2 0.4 0.6 0.8 1 

EOC(\\e\\ 2 o, m ) 



0.2 0.4 0.6 0.8 1 

IEI(E t 2 m + E 3 2 m) 



"i 1 1 r 



J I I L 



"i 1 1 r 



j i i i_ 



i 1 1 r 



J I L 



"i — i — i — r 



0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



0.2 0.4 0.6 0.8 1 



(b) Example (81) solution with fixed time step A ~ h . Inverse effectivity index is asymptotically tending towards value 
0.004. 



Figure 4: Example (81) with solution U2, L°° estimators and error. 




(b) Example (81) solution with fixed time step A w h • Inverse effectivity index is asymptotically tending towards a value 
around 0.005. Note the improvement on the IEI behaviour as the time step is refined. 



Figure 5: Example (80) with solution iai, adaptive algorithms based on L 2 estimators. 




Figure 6: Example (80) with solution m, adaptive algorithms based on L 00 estimators. 




(b) ImplicitTimeStepControl algorithm in comparison to the uniform approach converging to the L°° error tolerance 
pa 0.37. Depicted from the left are: timestep length, accumulated degrees of freedom and error evolution over time. 



Figure 7: Example (81) with solution U2, adaptive algorithms based on L 2 estimators. 




(b) ImplicitTimeStepControl algorithm in comparison to the uniform approach converging to the L 2 error tolerance 
0.0026. Depicted from the left are: timestep length, accumulated degrees of freedom and error evolution over time. 



Figure 8: Example (81) with solution 112, adaptive algorithms based on L°° estimators. 




(b) ImplicitTimeStepControl algorithm in comparison to the uniform approach converging to the L°° error tolerance 
0.0037. Depicted from the left are: timestep length, accumulated degrees of freedom and error evolution over time. 



